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The electronic spectra of rotationally faulted graphene bilayers are calculated using a continuum 
formulation for small fault angles that identifies two distinct electronic states of the coupled sys- 
tem. The low energy spectra of one state features a Fermi velocity reduction which ultimately 
leads to pairwise annihilation and regeneration of its low energy Dirac nodes. The physics in the 
complementary state is controlled by pseudospin selection rules that prevent a Fermi velocity renor- 
malization and produce second generation symmetry-protected Dirac singularities in the spectrum. 
These results are compared with previous theoretical analyses and with experimental data. 



PACS numbers: 73.22.Pr, 77.55. Fx, 73.20.-r 

The variation of the electronic properties of few layer 
graphenes (FLG's) with their layer stacking is receiving 
increasing attention. FLG's represent a family of materi- 
als that bridge the pseudo relativistic properties of single 
layer graphene with the more conventional semimetal- 
lic behavior of bulk graphite. The atomic registry of 
neighboring layers and stacking sequence are structural 
parameters that determine their electronic propertie^^^. 
In twisted FLG's where the crystallographic axes of 
neighboring layers are misaligned by a rotation angle 
titt/S the interlayer interactions produce remarkably 
rich physics that is being actively studiecpHU 

This paper presents a continuum theory of the low en- 
ergy electronic physics in twisted bilayer graphenes for 
small rotation angles, as illustrated in Fig. 1. Our ap- 
proach reveals the existence of two distinct electronic 
states in this system that present quite different elec- 
tronic properties. The behavior of one state is identi- 
fied with the situation described by a frequently adopted 
continuum formulation of this problem® the inter- 
layer coupling renormalizes the Fermi velocities of the 
individual layers and hybridizes their Dirac cones in the 
spectral region where they merge. In the complemen- 
tary state we find that the Fermi velocity renormaliza- 
tion is nearly completely prevented by a pseudospin se- 
lection rule and the interlayer hybridization inherits a 
novel momentum space geometry producing a set of sec- 
ond generation Dirac singularities. The behavior in this 
latter family agrees well with properties experimentally 
observed for rotationally faulted FLG's thermally grown 
on SiG (0001)^ suggesting that this physics is real- 
ized in this form of FLG. We briefly discuss the relation of 
our new results to prior theoretical and to experimental 
studies of these systems. 

The physics described below is identified by consider- 
ation of the effects of the lattice symmetry on the low 
energy electronic physics. We show that the geomet- 
rical structure of the low spectrum is determined by a 
symmetry- allowed threefold anisotropy in the interlayer 
coupling amplitudes which, though absent from conven- 
tional two-center tight binding models, occur in empirical 
models of interlayer coupling in graphite. We find that 




FIG. 1: Lattice structures of two twisted graphene bilayers 
rotated away from AA stacking by angles = 3.89° (top) 
and Q — 56.11° (bottom). The insets show schematically the 
dispersions of two nearby Dirac cones in these structures in 
the absence of their interlayer coupling. 



the sign of this anisotropy distinguishes two quite differ- 
ent electronic states of this system. 

The coupling between the two sublattices in the two 
layers can be represented by a (position dependent) 2x2 
matrix operator Ti2{r). As shown in Figure 2, for small 
angle faults the registry between layers in the unit cell 
evolves smoothly from regions locally characterized by 
AB (region a), BA (/3) and AA (7). The smoothest 
possible supercell-periodic matrix-valued expression for 
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Coefficient Parameterization I 



II 



tg 



Cab 



(71-73)79 43.3 8.3 
74 + (71 - 73)/3 130.0 69.0 
(71 + 273)73 130.0 340.0 



TABLE I: Fourier coefficients (meV units) for the interlayer 
hopping operator Eqn. 2, fitted to the Slonczewski Weiss 
McClure parameterization for Bernal stacked layers. Model 
I: 71 = 390 meV, 73 = 74 = 0. Model II: 71 = 390 meV, 73 
— 315 meV and 74 = 44 meV. 



T(f) is given by the expansion 

6 

with constant matrix coefRcients in and where Qn are the 
six elements of the first star of reciprocal lattice vectors 
dual to the superlattice translations Ti and T2 . The ma- 
trix coefficients in (n = 1,6) can be determined from the 
couplings in the locally registered regions; for example in 
the geometry of Figure 2 the even elements of the first 
star have coefficients 



in c 



= tg 



z 1 

z* z 



(2) 



where z = e^'^^l'^ and the coefficients for the odd elements 



are tn odd — 



The constant matrix has the form 



Cba Chb 



(3) 



with real coefficients satisfying Caa = Cbf, and Cab = 
Cba- The interlayer operator of Eqns. (2) and (3) is 
thus parameterized by three real constants tg, Caa and 
Cab- We choose these coefficients so that the interlayer 
matrix T{fa) matches the Slonczewski- Weiss-McClure 
(SWMcC) interlayer parameters 71, 73 and 74 for Bernal 
stacked graphite shown in the inset of Figure with the 
results in Table I. We note that the 73 parameter (hop- 
ping between unaligned sublattices in the two layers) is 
comparable to 71 and that the 74 parameter (hopping be- 
tween aligned and unaligned sublattice sites) is relatively 
weak. 

The conven tiona l continuum description of twisted bi- 
layer graphen^^E^ can be derived from the constant ma- 
trix tg . The low energy Hamiltonian is a long wavelength 
expansion around the zone corner points in each layer; in 
this Dirac basis the matrix elements in Eqn. 1 acquire 
the phases exp(i(G' ■ t[ — G ■ fj)) where G{G') are recip- 
rocal lattice vectors in the two separate layers and Tj{fl) 
are sublattice positions. Boosts by a triad of (G, G') 
pairs translate the Hamiltonian to three pairs of zone 
corner points that are separated by Ai? and its ±27r/3- 
rotated counterparts. This generates three possible con- 
stant coupling matrices indexed by the momentum dif- 
ferences AT?;. With a conventional choice of origirP^^ 




FIG. 2: Lattice structure for a segment of twisted bilayer at 
rotation angle 3.89°, with superlattice translation vectors T\ 
and T2. The points labelled a, /3 and 7 are high symmetry 
points in the unit cell. The inselP^ illustrates three hopping 
processes in the interlayer Hamiltonian. 



where the A sublattice site of one layer is aligned with 
the B sublattice of the other, the matrices are 



Ti = 




^aa ^ab \ 
Cba Ctb J 

^Caa Cab 
Z*Cba ZCbb 



z c„ 



Cab 



ZCba Z Cbb 



(4) 



In one of these valleys the Hamiltonian for the coupled 
bilayer with a momentum offset AK is 



jj ^ ( fivpcr ■ (-iV) 



Tl 

hvpCTe ■ (-iV - AK) 



(5) 



where ag are Pauli matrices resolved along the axes of the 
(^-rotated layer. The problem can be written dimension- 
less form by scaling all momenta by the offset \AK\ and 
energies by Eg — hvF\AK\. The scaled coupling coefh- 
cients are c = c/Eg = 3ac/{8T:hvF sm{6/2)) (where a is 
the single layer graphene lattice constant) which increase 
with decreasing rotation angle. 

Model I (Table I) is an isotropic interlayer model with 
73 = 74 = 0. For an isotropic coupling model Caa — Cbb = 
w and the interlayer matrices are 



1 1 
1 1 



■,T2 = w 



z 1 

z* z 



z* 1 

z z* 



(6) 
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FIG. 3: Electronic spectra for twisted bilayers using the 
interaction parameters (a): to — c(l + a^), c = 0.21, (b): 
to ~ cax, c — 0.55 and (c): to — cl,c = 0.55. Qpar and Qpcrp 
are momenta in units of the offset lAi^l and the ordinate is 
the scaled energy E/Eg — E / {hvpAK). In (b) Dirac cones 
with opposite helicity are coupled, in (c) Dirac cones with 
the same helicity are coupled. The insets give the locations 
of singular points in the spectra describing the annihilation 
and regeneration of Dirac nodes (red diamond) in the com- 
pensated case (b) and the appearance of a singular point of 
degeneracy (C) for the uncompensated case (c). The point C 
is a second generation Dirac point singularity in the coupled 
spectrum. 



with w = 130 meV. The form of these matrices and their 
prefactor agree with the estimates ( w ~ 110 meV) ob- 
tained from tight binding calculation^^EH, Our construc- 
tion shows that these terms project the q = term of the 
interlayer potential into the Dirac -fC-point (pseudospin) 
basis thereby coupling the electronic states in the two 
layers with identical crystal momenta. Since only the 
q = term in the coupling is retained it does not de- 
pend on a relative lateral translation of the two layers, 
in agreement with earher work^^ and physically reason- 
able since for small twist angle a rigid layer translation 
produces insignificant changes to the Moire superlattice. 
Thus Model I reproduces the existing continuum theo- 
retic phenomenology of the coupled system, and the cal- 
culation leading to Eqn. 6 provides an alternate (and 
compact) derivatio n of t he effective Hamiltonian used in 
these earlier studiepl^. The left panel of Fig. 3 shows 
the bilayer spectra computed in this model which shows 
the expected (^-dependent) reduction of the Dirac cone 
velocities and a hybridization of the two branches in the 
spectral region where they merge. 

We now consider a refinement of the interlayer cou- 
pling matrices using the parameterization of Model II. 
The salient properties of the SWMcC parameterization 
are the introduction of the interlayer amplitudes 73 and 



74 with 73 comparable to 71 and 74 significantly smaller. 
Note that 73 and 74 represent interlayer hopping pro- 
cesses at the same range but in different directions with 
respect to the layer crystallographic axes. The asym- 
metry between 73 and 74 thus reflects an intrinsic three- 
fold lattice anisotropy in the interlayer amplitudes which, 
though symmetry-allowed, does not occur in the isotropic 
two center tight binding approximation. Significantly, 
these additional terms break the symmetry between the 
pseudospin-diagonal and off diagonal terms in to (Table 
I) so that tlie coupling matrix is dominated by its off di- 
agonal amplitudes. An instructive limit considers to cx ax 
for which the Fig. 3(b) shows the spectrum calculated 
for a 6* = 3.89° rotation away from Bernal stacking. Here 
the two Dirac cones have merged at low energy producing 
two composite low energy singular points. Note that the 
linear low energy dispersion is replaced by an approxi- 
mately quadratic form near the center of symmetry of 
these spectra and that the momentum offset between the 
singular points in the spectrum is along the Qporp axis, 
i.e. 7r/2-rotated with respect to the original Dirac cone 
offset AK. 

These spectral changes reflect the proximity to a criti- 
cal point that occurs at c = 1/2 in this theory. This can 
be understood by considering a single layer sublattice 
exchange operation implemented by the gauge transfor- 
mation 

^ ^ / I 0\( HK{q) CO. W I \ 

\Qox)\ ~cGx HK{q- AK) J \ a, ) 





cl 


cl 


■ Hxiq- AK) ■ &x 



demonstrating that this system has a scalar coupling 
Dirac cones with compensating helicities (Berry's phase 
±7r). Increasing the control parameter c (by decreasing 
6) draws the nodes together until they become coincident 
at a critical coupling strength c = 1/2 and annihilate 
(Fig. 3b inset). For c > 1/2 new singularities emerge 
at £^ = separated by AQ directed perpendicular to the 
original offset AK. Using the parameters listed in Table 
II, c{9 — 3.89°) — 0.55, i.e. just on the strong coupling 
side of this transition. The residual curvature in the low 
energy spectrum and the associated reorientation of AQ 
are both clearly evident in Fig. 3b. It is noteworthy 
that the momentum separation between the zero energy 
contact points is generally not determined purely geomet- 
rically by the rotation angle as is generally assumed, but 
instead is modified by the interlayer coupling. This oc- 
curs because the interactions between layers produces an 
effective gauge field seen within each layer that shifts the 
momentum of its zero energy states. The tt/2 rotation 
of the momentum offset that bridges the contact points 
on the strong coupling side of the transition is a striking 
consequence of this gauge coupling. 

Reversing the sign of the threefold anisotropy in the 
interlayer matrix to produces a distinct electronic state. 
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The complementary behavior is understood by consid- 
ering the limit to on I which describes the coupling of 
Dirac cones with the same helicity, preventing annihi- 
lation of the Dirac nodes and leading to a qualitatively 
different geometry in the bilayer spectrum (Fig. 3c). The 
dispersing bands from the uncoupled cones are degener- 
ate everywhere along the line that bisects AK. How- 
ever, along the line that connects the Dirac nodes the 
pseudospins are orthogonal and the interlayer coupling 
is symmetry-forbidden, turning on linearly as a function 
of the transverse momentum Qperp- Thus the coupled 
system retains a twofold point degeneracy midway be- 
tween the displaced Dirac noded^SJ ^jjg cancellation of 
the interlayer coupling at this critical point is the bilayer 
analog of the "absence of backscattering" due to the tt 
Berry's phase in single layer graphene. In the vicinity of 
this critical point interlayer coupling is allowed and pro- 
portional to the transverse momentum. Thus this system 
exhibits second generation Dirac singularities in its cou- 
pled layer spectrum as shown in Fig. 3c: hybridization of 
the two layers is symmetry forbidden at a discrete critical 
crossing point. We refer to this complementary state as 
the uncompensated bilayer state. 

The relative helicity of the two Dirac cones also con- 
trols the renormalization of their Fermi velocities, further 
distinguishing these two states. For Dirac cones of oppo- 
site helicities, perturbation theory on the Hamiltonian in 
Eqn. 7 for small c modifies the velocity operators 

v+ — vpa+ vp{l — (?)(7+ 

i)_ — vf<7- — > wf(1 — c^)a- (8) 

which symmetrically reduces both and Vy ; summation 
over the triad of offset momenta AKi yields the renormal- 
ized velocity Vp = vp{l — 9c?) exactly as found in earlier 
work^ By contrast for coupling between nearby cones 
of the same helicity perturbation theory yields 

■0+ = vpd^ — > vp{(T^ — <?a^) 

u_ — vpa^ vp{a^ — (?'y+) (9) 

so that in this class the corrections to the velocity are 
weaker, oc c^. Moreover they have a twofold cos(20) 
anisotropy so they vanish by symmetry after summing 
over the threefold symmetric triad of AKi. Thus the 
Fermi velocity is unchanged by the interlayer coupling in 
this class of bilayers. Physically this can be understood 
by observing that the bands dispersing through the Dirac 
nodes are connected smoothly to the second generation 
points of degneracy at AK /2. 

The distinction between the compensated and uncom- 
pensated states in the small angle limit reflects a lattice- 
scale property that determines the matrix structure of 
the long wavelength coupling in Eqn. 1. This should 
be distinguished from the different mechanism by which 
sublattice exchange symmetry determines the direct cou- 
pling between the Dirac node^I^. The latter requires 
finite momentum umklapp interlayer hopping processes 



which, though significant for low order rational commen- 
surate rotations, are negligible in the small angle limit 
considered here. For example, note that sublattice ex- 
change "even" and "odd" commensurations are related 
by a rigid sublattice translation of one layer at a fixed 
rotation angle. In the small angle regime this translation 
simply permutes regions of the bilayer that are locally 
in AB, BA and A A registry as shown in Fig. 1, but 
it does not change tg which determines the spectrum. 
Thus sublattice exchange "even" and "odd" structures 
become indistinguishable in the small angle limit. Note 
also that bilayers at rotation angles 6 and 9 — tt/3 — 6 are 
commensuration pairs that can be distinguished by their 
sublattice exchange paritj^^. Even and odd parity com- 
mensurations are, respectively, inflated generalizations of 
the primitive AA and AB stacked bilayers. This symme- 
try ultimately determines the valley structure of the in- 
terlayer amplitudes that directly couple the Dirac nodes 
of neighboring layers. This interlayer umklapp coupling 
derives from the finite momentum terms in the interlayer 
Hamilonian in contrast to the q = terms that control 
the physics for small angle rotations. 

The spectra for these two classes are ultimately deter- 
mined by the pseudospin asymmetry in t^. The conven- 
tional SWMcC model selects the class that couples cones 
with compensating helicities. In this model the spec- 
tral transition illustrated in Figure 3 occurs for rotation 
angles near 4°, i. e. in a range that is frequently stud- 
ied experimentalljEEH xhe physics of the uncompen- 
sated class occurs for Caa > Cab which requires 74 > 73. 
Although this is excluded by the conventional SWMcC 
parameterization it is important to note that this pa- 
rameterization is designed to fit data for Bernal stack- 
ing, and it likely does not properly represent the matrix 
structure of the coupling in AA registered regions. In 
particular using the parameterization of Table I, the spa- 
tial dependence of Eqn. 1 shows that strong interlayer 
coupling in A A stacked regions requires 74 > 73. Micro- 
scopically this originates from interlayer tunneling pro- 
cesses along the edges of eclipsed hexagons in the aligned 
AA structure, a motif which does not appear at all for 
Bernal stacking. In the spirit of the SWMcC theory it 
is therefore appropriate to retain 73 and 74 as param- 
eters which can be determined from the experimentally 
observed properties of twisted graphenes. 

In fact the phenomenology of the uncompensated class 
provides a striking explanation for many of the puz- 
zling observed spectral properties for rotat ionally faulted 
graphenes thermally grown on SiC(000l )^^l^^l^^ [ Lan- 
dau level spectroscopy shows a negligible renormaliza- 
tion of the Fermi velocity in these structureJi^ and fur- 
thermore angle resolved photoemission finds no evidence 
for a hybridization-induced avoided crossing of the inter- 
secting Dirac cones, despite a careful searclP^. This is 
completely consistent with the existence of a node in the 
interlayer coupling at the midpoint between offset Dirac 
cones characteristic of the uncompensated class. This as- 
signment can be confirmed definitively by measurements 
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of the quasiparticle dispersion along an azimuth passing 
through the midpoint between the displaced Dirac cones, 
but perpendicular to AK; these should show a band split- 
ting linear in the transverse momentum around the point 
of degeneracy. Alternatively, if these bilayers exist in the 
compensated class, photoemission should be able to de- 
tect the annihilation and re-emergence of their singular 
contact points along with the band curvature in their 
spectra in the crossover regime as illustrated in Fig. 3(b). 

By contrast, experiments on rotationally faulted CVD- 
grown graphenes have observed phenomena that have 
been associat ed wit h the spectral properties of the com- 
pensated clas d^^ l ^° [ Features due to the van Hove sin- 
gularities arising from an the avoided crossing of hy- 
bridized Dirac conefJIS' and a ^-dependent low energy ve- 
locity renormalization have both been reported^^". These 
features are at least qualitatively consistent with the pre- 
dicted behavior of the compensated class and have been 
analyzed within a theoretical model representative of this 
clasa^. We note that these measurements study sam- 
ples at small rotation angle where the proximity to the 
merger of the Dirac singularities (Fig. 3) should be man- 
ifest in these data though their effects have not yet been 
considered in the analysis. It is interesting that these 
samples exhibit a large periodic height modulation w lA 
in the superlattice unit cell peaked in the AA-registered 
zones'^® . It is tempting to speculate that these CVD sam- 
ples are grown as rippled structures that partially delami- 
nate in these regions thereby locally weakening their con- 
tribution to the q = coupling coefficients. In this sce- 



nario the strongly coupled regions would maintain Bernal 
registry as described by the conventional SWMcC param- 
eterization and identify these samples as members of the 
compensated bilayer family. 

The distinction between the two complementary states 
is controlled by an important three-fold anisotropy in 
the interlayer tunneling amplitudes. This physics is not 
captured by an isotropic two-center tight binding theory, 
which inevitably leads one to the coupling model in Eqn. 
6 which happens to occur at a crossover between two 
rather different electronic models for the system. The 
effects of the threefold anisotropy are accessible in den- 
sity functional calculations of these structures, but for 
practical reasons these have been restricted to short pe- 
riod superlattices which do not address the small angle 
regime where the continuum theory is most appropriate. 
For short period commensurate structures, the Fermi ve- 
locities found in these calculations are consistent with 
the values for single layer graphene. This could arise 
from the small value of c in the large angle regime, the 
intrinsic behavior of the uncompensated class or an in- 
terlayer mass term which is important for short period 
superlattices^''. 
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